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We show that an equilibrium-like additivity property can remarkably lead to power-law distribu¬ 
tions observed frequently in a wide class of out-of-equilibrium systems. The additivity property 
can determine the full scaling form of the distribution functions and the associated exponents. 

The asymptotic behaviour of these distributions is solely governed by branch-cut singularity in the 
variance of subsystem mass. To substantiate these claims, we explicitly calculate, using the addi¬ 
tivity property, subsystem mass distributions in a wide class of previously studied mass aggregation 
models as well as in their variants. These results could help in thermodynamic characterization of 
nonequilibrium critical phenomena. 
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I. INTRODUCTION 


Simple power-law scaling is ubiquitous in nature [ij. 
They ^pear in the distributions of drainage area of 
rivers [3, droplet size [^, 3, size of clusters formed in 
polymerization processes rain size size of fr^- 
ments in fractured solids [? , population and wealth [1,[9[, 
and in stock market fluctuations 0, etc. Evidently, 
power laws, which are usually associated with criticality 
through emergence of a diverging length scale, are ob¬ 
served in widely unrelated systems, suggesting existence 
of some broad underlying principle. Recent evidence that 
living systems might be operating, independent of most 
of the microscopic details, in the vicinity of a critical 
regime m indeed invoke further questions - how and 
why systems adapt to near-criticality. 

There have been several attempts to reveal the origin of 
the power laws in nature, through studies of paradigmatic 
nonequilibrium models - most appealing being sandpile 
[12 - [i 3| and mass aggregation models [isl - l^ . Many of 
these models - where there is a conservation law or, in 
case of violation, the law is weakly violated in the sense 
that the systems are slowly driven - are intimately con¬ 
nected to each other. For example, the mass aggrega¬ 
tion models [IMlllIl are connected to directed abelian 
sandpile model 221 or to the models of river network Q. 

In this paper, we argue that power-law distributions 
in out-of-equilibrium systems can arise simply from ad¬ 
ditivity property, the tenet of equilibrium thermodynam¬ 
ics. We find that the divergence in the response function 
is the key: Diverging fluctuations can, in principle, arise 
from distributions other than power laws, which are how¬ 
ever prohibited if one imposes additivity and consequent 
fluctuation-response (FR) relation. The response func¬ 
tion determines the full scaling form of the distribution, 
at as well as away from criticality, and critical exponents 
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originate from the singularity in the response function. 
To demonstrate this, we consider mass aggregation mod¬ 
els which are known to have nonequilibrium steady state 
with scale invariant structures. At all mass densities, the 
distribution function (m) of mass m in a subsystem of 
volume V, which is obtained solely from the FR relation, 
is shown to have a scaling form Py{m) ~ m~'^ exp(flm). 
The quantity p,{p) = p,{p) — p.{pc), inverse of a cut-off 
mass m*{p) = — l//2(p), is an analogue of equilibrium-like 
chemical potential and provides a useful thermodynamic 
interpretation of the emergence of power laws in nonequi¬ 
librium steady states. The exponent r and the critical 
properties of chemical p{p) arise from a multiple-pole or 
branch-cut singularity in the variance at a critical mass 
density pc- As the critical density is approached p ^ pc, 
nonequilibrium chemical potential vanishes p,{p) —> 0, 
leading to pure power laws. Beyond the critical density 
p > Pc, there is a gas-liquid like phase coexistence. 

The above result immediately provides answer to why 
power law, at or away from criticality, appears 
so often in mass aggregation models - especially in 
higher dimensions, at all densities and irrespective of 
that the motion of the diffusing masses is biased or not 
[H)[ll|,[2l,[23|- Interestingly, the same power law appears 
in fc-mer distribution in the classic Flory-Stockmayer 
theory of polymerization and also in particle number dis¬ 
tribution in three dimensional ideal Bose gas near crit¬ 
ical point, irrespective of whether the systems are in or 
out of equilibrium - thus indicating a universality. We 
demonstrate that the law is a consequence of a 

simple-pole singularity in the variance. The whole anal¬ 
ysis is extended also to nonconserved mass aggregation 
models. We validate our theory by explicitly calculating 
mass distributions in previously studied mass aggrega¬ 
tion models and their variants and by comparing them 
with simulations. 

Organization of the paper is as follows. In section II. A, 
we discuss additivity property; in section II.B, we discuss 
the connection between singularity in the variance and 
the asymptotic behaviour of the mass distribution func- 
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tion. In section III, we illustrate our analytic methods 
in a broad class of model systems, both in conserved- 
mass aggregation models and its nonconserved versions. 
Finally, we summarize our results with a concluding per¬ 
spective. 


II. THEORY 
A. Additivity Property 

We start by invoking an additivity property which a 
wide class of systems, irrespective of that they are in or 
out of equilibrium, could possess. Consider a continuous- 
mass (generalization to discrete masses is straightfor¬ 
ward) transport process on a lattice of size V and then 
divide the system in i/ = V/v identical subsystems or 
cells, kth subsystem with mass m^, where total mass 
M = remains conserved. Provided the subsys¬ 

tems are large compared to spatial correlation length, 
they could be considered statistically almost independent 
[27H30j| . In that case, the joint subsystem mass distribu¬ 
tion in the steady state can be written in a product form. 


MS’"*-" 


Z{M, V) 


( 1 ) 


where Wy(rnk) is an unknown weight factor (to be deter¬ 
mined later; see Eq. m depending only on the subsys¬ 
tem mass ruk, 2' = H/c / dmkWy{mk)6{J2kknk - M) = 
exp[—Vf{p)] the partition sum, f{p) a nonequilibrium 
free energy density and p = M/V mass density. The 
product form in Eq. [T] amounts to an equilibrium-like 
additivity property, in the sense that a free energy func¬ 
tion F = '^f, in Wv{mk) is minimized in the macrostate. 

Using standard statistical mechanics [1^, Eq. [T] leads 
to the probability distribution Prob[mfc € (m, m+dm)] = 
Py(m)dm for subsystem mass where 


Py{m) = 


,{m)t 


( 2 ) 


with p,{p) a nonequilibrium chemical potential, and Z 
the normalization constant. The weight factor Wy{m) 
and chemical potential p{p) = df /dp can be obtained 
using a fluctuation-response relation [27h31I |. 


dp 

dp 




( 3 ) 


where the scaled variance cr^(p) = (I/u)((m|) — (nik)'^) 
in the limit of 1. The free energy density function 
f{p) can be obtained through the relation p{p) = df /dp, 
i.e., f{p) = f p{p)dp + /3 with chemical potential p{p) = 
f l/cr^(p)dp + a (obtained from Eq. [3]) and a and /3 arbi¬ 
trary constants of integration. Then Laplace transform 
of Wy{m) is written as Wy{s) = f Wy{m) exp{—sm)dm = 


i.e.. 


„-A„(s) _ 


( 4 ) 


Then, the function Xy{s) can be obtained from Leg endre 
transform of free energy density function /(p) [32|, 

Xy{s) = V [infp{/(p) -h sp}] = v[f{p*) + sp*], (5) 

where p*(s) is the solution of 


s = -p(p*)- 


( 6 ) 


As discussed later, Eq. [5] requires concavity and differen¬ 
tiability of /(p). In the discrete case, the weight factor 
can be calculated as Wy/m) = (I/27ri) Wy(z)/z^~^^dz 

where Wy{z) = J2m=o F^Wy/m) is obtained from w[s) by 
substituting s = — In z with C a suitably chosen contour 
in the complex z-plane. 


B. Singularity in Variance and Mass Distribution 

Importantly, the distribution function Py(m) is deter¬ 
mined solely by the functional form of the scaled variance 
tT^(p). We argue below that singular response functions 
generate only power-law distributions. Other functional 
form of mass distribution Py (rn) with diverging variance 
is also possible [s^, which, we show, however are not al¬ 
lowed if the FR relation holds. In this paper, we mainly 
focus on multi-pole singularity at a finite density pc, 


I sip) 

a^/p) = (Pc-p)" 

oo 


for p < Pc, 
otherwise. 


( 7 ) 


This form, with 0 < n < oo, is relevant in the con¬ 
text of a wide class of mass aggregation models as dis¬ 
cussed in section III. The analytic part p(p) is not par¬ 
ticularly relevant in determining the asymptotic form of 
the distribution Py{m), however it contributes to the ex¬ 
act form of Py{m) (discussed in section II.B.IV). In fact, 
other kinds of singularities, such as logarithmic singu¬ 
larity cr^(p) ~ [ln(pc — p)]^ or exponential singularity 
exp[(pc — p)~s] where p > 0, I/|p — pd”, and the case 
with n < 0 can also arise. One can show that they all 
lead to power laws, possibly with logarithmic corrections 
to the power-law scaling (discussed in the following sec¬ 
tions). 

The divergence in the variance, as in Eq. [7] (or in 
the cases of logarithmic and exponential divergence), in¬ 
deed has broad implications, not only in conserved-mass 
aggregation models but also in their nonconserved ver¬ 
sions. Note that the FR relation in Eq. implies that 
free energy density /(p) is not a strictly concave function 
of p and has a linear branch of slope p(pc) for p > Pc- 
Moreover, f”{p = pc) = p’[p = pf) = 0 (prime denotes 
derivative w.r.t. p) implies a point of inflection in / — p 
curve at p = Pc- That is, free energy density function 
can be written as 


fix) = 


Wy{m)e ^"^dm. 


f p(x)dx for X < Pc 

piPc)ix - Pc) + /(Pc) otherwise. 


( 8 ) 
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FIG. 1: (Color online) Schematic representation of conden¬ 
sation transition. Panel (a): variance o-^{p) as a function of 
density p; panel (b): chemical potential /r(p) as a function 
of p; panel (c): free energy density /(p) as a function of p; 
panel (d): Legendre transform A(s) of free energy density as 
a function of s. 


Q,)i+i/(i+")] in leading orders of (s -I- a), implying 


g-am 

Wy{m) -—, 

m/ 

where, for large subsystem masses m 3> u, the power-law 
exponent t in the denominator is given by 



1 

(TT^J ’ 


( 12 ) 


with the following inequality 2 < r < 3 (since 0 < n < oo 
in Eq. ED. This translates into the mass distribution 
having a scaling form, 


Pv{m) oc 


mJ 




(13) 


where jl{p) = l/tT^(p)dp = ^(p) — p.{pc) is an effec¬ 
tive chemical potential, inverse of which gives a cut-off 
m* = — l/p in the distribution, and the scaling func¬ 
tion <i>(a::) = exp(—a;). Later, we explicitly calculate 
p{p) in specific model systems. Note that /i(pc) = 0 at 
critical point p = p^ and consequently Py(m) becomes 
a pure power law. Moreover, by defining a critical ex¬ 
ponent 6 = 1 -I- n as p{p) ^ {pc — pY, we get a scaling 
relation S{t — 2) = 1. 


Consequently, Legendre transform of f{p) develops a 
branch-cut singularity (see Eq. 1111) : for schematic repre¬ 
sentation of the above analysis, see Fig. El This construc¬ 
tion of a nonequilibrium free energy function f{p) from 
a general thermodynamic consideration readily explains 
the phase coexistence between a fluid and a condensate, 
as observed in the past in many out-of-equilibrium sys¬ 
tems (discussed in section III). 


2. Logarithmic singularity 

Now, we consider the case of logarithmic singularity 
where variance cf‘^{p) diverges logarithmically as given 
below, 

0-2(p) = / 9{p)Y^{Pc - P)Y for p < Pc, 

y oo otherwise. ' ' 


1. Multi-pole singularity 


Integrating Eq. [3] near p = pc, we obtain chemical po¬ 
tential, in leading order of (pc — p). 


To analyse the behaviour of A„(s) in the case of Eq. El 
we integrate Eq. [3] near p = Pc and obtain 


p(p) ~ 


{Pc - p) 

g{pc)[\n{pc - p)]P 


H” Q:, 


(15) 


(f) _ 

liw 1 (9) 

{n + l)g{pc) 

which gives (pc — p*) — [(n -I- l)p(pc)(s + by 

using Eqs. IHandlE] Integrating chemical potential p(p), 
we get free energy function 


f{p) 


jPc - p)"+^ 

(n -h l){n 2)g{pc) 


ap-\- (3 


( 10 ) 


and write Xv{s) = v[f{p*) -I- sp*], in leading order, as 


Xy{s) ~ v[ao -I- ai(s -I- a) + 02(3 -I- q;)"+i]. 


( 11 ) 


which gives (pc — p*) — < 7 (Pc)[ln(pc —p*)]^(s + Q:) from Eq. 
El Free energy density is obtained by integrating above 
chemical potential p(p). 


f{p) - 


{Pc - pf 

2p(pc)[ln(pc - p)]P 


ap-\- j3, 


(16) 


and, accordingly, its Legendre transform, in leading or¬ 
ders. 


Xy{s) ~ u[ao -l-ai(s-|- a) + 02 ( 3 -!- a)^{ln(s-|- a)}^]. (17) 


For large mass v, this implies that the weight factor 
has a functional form of a power law with logarithmic 
correction. 


where Qq, 01,02 are constants. Thus, we obtain w{s) = 
exp[—A„(s)] ~ const, x [1 — vai{s -I- a) — va 2 {s -\- 


Wv{m) 


(Inm)^’ ^ 
rrP 
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and the corresponding mass distribution function, 


Pvijn) oc 


(In m)P ^ 
w? 




(18) 


where effective chemical potential jl{p) = p,{p) — p{pc) 
and the power-law exponent in the denominator is r = 3, 
the borderline case of Eqs. 171 and IT^ with n = 0. 


where the denominator is essentially a sub-leading cor¬ 
rection to the free energy function, with a ^ 0(v) being 
a model-dependent cut-off mass. The correction term is 
obtained from the fact that, for large masses m ^ v, free 
energy function f{x) has a linear branch (see Eq. [8]) and 
the mass distribution function must have the asymptotic 
form Py{m) ~ exp(flm) as in Eq. 1131 


3. Exponential singularity 


III. MODELS AND ILLUSTRATION 


We also consider the case where the variance diverges 
exponentially, a^{p) ^ exp(pc — p)~^ for P < Pc and 
(j‘^(p) = oo otherwise. The analysis is similar to the 
ones given above. Substituting chemical potential p{p) ~ 
const X (pc — pY^^ exp [—(pc — p)~Y + a in Eq. [5] and 
solving in leading order of (s -I- a), we get (pc — p*) ^ 
{ln(s -b and consequently 


\y{s)c:iv Oo-b ai(s-b a)-b 02 ( 5 -b a){ln(s-b a)} 


(19) 

For large mass v, this leads to the mass distribution 
function 


(Inm) ^ 
Pv{m) oc 1 




( 20 ) 


where p,{p) = p(p) — p(pc) and the power-law exponent 
in the denominator is r = 2, the borderline case of Eqs. 
[7l and [T2l with n = 00 . 


4- Subsystem mass distribution 

For any finite v, it is not easy to find the distribution 
function Py (m) at small or intermediate values of masses 
m V because, in that case, one requires to invert Eq. 
m using inverse Laplace transform, i.e., by evaluating the 
integral. 


Wv{m) = ^. j ( 21 ) 

27rz Jc 

on the complex s-plane; the contour C on the complex 
plane should be chosen such that the integral converges. 
However, in the models we consider here, it is not possible 
to get an exact functional form of A^(s) for all s, which 
actually involves solving the transcendental Eq. [51 How¬ 
ever, for large subsystem size u ^ 1, the analysis sim- 
plies as the function —{1/v) In ^^(to) is simply related to 
Xy{s)/v by Legendre transformation [s^ and therefore, 
in leading order of m, is the free energy density function 
f{m/v) itself, which has been already constructed in Eq. 
[5] (see Fig. [T] and related discussions). The subsystem 
mass distribution function now can be written as 


Pv{m) cx 


^-vf(m/v)+fi{p)m 


( 22 ) 


A. Conserved Mass Aggregation Models (CMAM) 

We now substantiate the above claims in a broad class 
of nonequilibrium models which were studied intensively 
in the last couple of decades. In this paper, we mainly 
focus on the models having multi-pole singularity in the 
variance. Other singularities, e.g., exponential or log¬ 
arithmic, are possible, however not as common as the 
multi-pole one. For the purpose of illustrations, we 
specifically consider the case with n = 1 and mass distri¬ 
bution at a single site, i.e., v = 1. We define a generalized 
version of conserved mass aggregation models (CMAM) 
studied in (l9l - [^ 1^ . for simplicity on a one dimen¬ 
sional lattice of L sites. We mainly focus on symmetric 
mass transfer, i.e., masses are transferred symmetrically 
to either of the neighbours. Here masses (or particles) 
diffuse, fragment and coalesce stochastically with either 
of the nearest-neighbour masses according to the follow¬ 
ing dynamical rules: (1) diffusion of mass rrn at site i to 
iztl with mass-dependent rate D{mi) where —>■ 0 and 
mi±i —>■ mi±i + rrii and (2) fragmentation of a discrete 
mass A at site z, provided A < m-i, and coalescence of 
the mass to either of the sites z ± 1 with mass-dependent 
rate w{A) where rrii —> to, — A and TOi±i —>■ -b A 

with A = 1, 2,... (continuous A considered later). Total 
mass M = conserved in this process. Even 

for these simple dynamical rules, the steady state weight 
in general is not exactly known. However, as spatial cor¬ 
relations are small, the additivity property as in Eq. |TJ 
to a good approximation, is expected to hold. 

We calculate the variance cr^(p) of mass at a single 
site in various special cases, using the additivity property 
Eq. [1] with V = 1. We take diffusion rate D{mi) = 1, 
independent of mass mj, w{A = 1) = wi (rate of single 
particle chipping), w{A = rrii — 1) = W 2 (rate of all-but- 
one particle chipping) and w{A) = 0 otherwise. 


1. Case I: CMAM with wi = 1,W2 = 0 

For wi = 1 and W 2 = 0 and symmetric mass transfer, 
the model becomes the symmetric one studied in ll^ 
(our model is a variant of those studied in [Ia[3)- For 
P A Pc, using additivity property, we exactly calculated 
the variance and consequently chemical potential with 
the critical density pc = \f2 — 1 We can calculate the 


(a -b mY'^ 
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FIG. 2: (Color online) Single-site (v = 1) mass distribution 
functions Pi(m) (points - simulations) in conserved-mass ag¬ 
gregation models (CMAMs) is compared with analytic expres¬ 
sion in Eq. [13] (lines - theory) for various densities. Panel (a): 
mass chipping rates wi = 1 ,W 2 = 0; panel (b): mass chipping 
rates wi = 0, UI 2 = 1. In all cases, mass diffusion rate D = 1 
and system size L = 5000. 


variance as given below (for details, see Appendix A) 

_ P(l + P)(l+P^) 

^ (/^/ /-I o 2 \ 

p{l + p){l + p^) 

[pc — p){V2 -I- 1 -I- p) 

with Pc = {\/2 — 1), for which one can obtain a chemi¬ 
cal potential p{p) and free energy function /(p), by in¬ 
tegrating the fluctuation-response relation as in Eq.2 in 
the main text. 


= J V + ln +a(24) 

and, upon one more integration. 


f{p) = J P{p)dp 


—2ptan ^p-bpln 




+ ap + (3 (25) 


calculated to be Pi{m) oc exp[(p(p) — p(pc))m]. 

(for details, see Appendix B). 

In panel (a) of Fig. |2l we have plotted single-site 
(u = 1) mass distribution function Pi{m), obtained from 
simulations, as a function of mass m for various values 
of densities p = 0.1, 0.2, 0.3, 0.414 and 1.0 and compare 
them with the theoretical expression in Eq. [T^l The 
theoretical results have a quite good agreement with the 
simulation results, especially at large mass to ^ 1. In 
panel (a) of Fig. O we have plotted subsystem mass 
distribution function P„(to), with v = 100, for densities 
p = 0.1 and 0.2 and compare them with the theoretical 
expression in Eq. I22l with cut-off mass a ~ 20; agreement 
between simulations and theory is reasonably good. 


2. Case II: CMAM with wi = 0,W2 = 1 


The CMAM with wi = 0, W 2 = 1 is a variant of the 
models studied in [ 2 ^. In this case, for p < Pc, the 
variance and chemical potential can be exactly obtained 
using additivity property (for details, see Appendix A). 
The variance is given by 

2 ^ ^ P(1 - P)(V - 2p+1) 

' = - 2f--4p+l - 

p(l-p)(2p^-2p-H) 

{pc- p){2 + V2-2p)' 

There is a simple pole at the critical density pc = 1 — 
1 / v^. By integrating fluctuation-response relation Eq. 
El we get chemical potential 


pip) 


2 tan ^ (1 — 2p) — In 


_2p(l - p) 


a, (27) 


and then free energy density 


/(p) = 2ptan ^(1 — 2p) — ln(l — p) + ln(l — 2p -|- 2p^) 


—pin 


1 


_2p(l - p) 


- 1 


ap + P (28) 


In panel (b) of Fig. [S] we have plotted single-site {v = 1) 
mass distribution function Pi (to), obtained from sim¬ 
ulations, for various values of densities p = 0.1, 0.15, 
0.2, 0.292 and 1.0. We find that the simulation results 
agree remarkably well with the analytical scaling form 
Pi (to) oc to“®/^ exp[(p(p) — p(pc))to] as in Eq. 1131 with 
T = 5/2 (for details of the derivation, see Appendix B). In 
panel (b) of Fig. jS] we have plotted subsystem mass dis¬ 
tribution function P„(to) for v = 100 for densities p = 0.1 
and 0.15 and compared them with theory, Eq. 1221 with 
a ~ 25; agreement between simulations and theory is 
reasonably good. 


3. Other Cases 


where a and P are two arbitrary constants of integration. We have also studied, through simulations, various 
For large mass m ^ 1, the mass distribution function is other cases (with P = 1): Case III. - w{A = I) = wi, 
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FIG. 3: (Color online) Subsystem mass distribution functions 
Pv{m) (points - simulations) in conserved mass aggregation 
models (CMAMs) is compared with analytic expression in 
Eq. [23 (lines - theory) for various densities. Panel (a): mass 
chipping rates uii = 1, u ;2 = 0 and a « 20; panel (b): mass 
chipping rates wi = 0, W 2 = 1 and cut-off mass a « 25. In 
all cases, mass diffusion rate D = 1, system size L = 10® and 
subsystem size v = 100. 



FIG. 4: (Color online) Variance vs. {pc — p). Black 

line is const, x {pc — p)~" with n — 1. Red rectangles are for 
Gase III (one and two particle fragmentation), blue triangles 
are for Case IV [A = 1,2,... is discrete with fragmentation 
rate w{A) = exp(—A)] and red circles are for case V [A > 
0 is continuous with fragmentation rate w(A) = exp( —A)]. 
Diffusion rate D{m) = 1 throughout. 


w{A = 2) = W 2 and w{A) = 0 otherwise. Case IV. - 
a discrete-mass model with w(A) = exp(—A) and Case 
V. - a continuous-mass model with w{A) = exp(—A). 
In these cases, in the absence of an analytical expres¬ 
sion of cr^(p), we checked in simulations (see Fig. |4]) that 
the variance near critical point indeed has the behaviour 
^ {pc — p)~‘^, with n = 1, which therefore leads to the 
same power-law exponent r = 5/2. 

One can also define an asymmetric version of the 
CMAMs discussed above. In one dimension, there are 
some nontrivial spatial correlations and the above mean- 
held analysis fails to capture the mass huctuations in the 
system. However, in higher dimensions, the above results 
qualitatively remain same also for the asymmetric mass 
transfer and is consistent with [ 2 ^. 

Interestingly, the exponent r = 5/2 appears also in 
the distribution of particle numbers in ideal Bose gas 
in three dimensions (3D) near the critical point where 
Bose-Einstein condensation (BEC) occurs. This could 
be easily understood from the fact that particle-number 
huctuation in the case of 3D Bose gas has the same crit¬ 
ical behaviour cr^(p) ~ {pc — with u = I, as in 

these ‘mean-held’ nonequilibrium systems having negli¬ 
gible spatial correlations. That, on a mean-held level, the 
nonequilibrium aggregation models belong to the univer¬ 
sality class of equilibrium Bose gas in 3D, so far has not 
been realized. 

It is quite instructive to consider a limiting case of 
Eq. [3 with n = 0, Pc = 00 and g{p) ^ p^~^ at large 
density, i.e., the variance cr^(p) ~ P^~^, with 6 < —1, 
diverges algebraically at inhnite density. As there is 
no singularity in the variance at any hnite density, our 
analysis quite straightforwardly shows that condensation 
transition cannot occur, consistent with the past observa¬ 
tions in the CMAM with mass-dependent diffusion [s^. 
Asymptotic scaling of the mass distribution can be ob¬ 
tained as follows. Using Eq. |31 we get p{p) ^ p^ (setting 
a = 0 without loss of any generality) and, consequently, 
Laplace transform of weight factor Wy{m), 

w„(s) ~ Oo + (29) 


immediately leading to mass distributions having a scal¬ 
ing form Py{m) oc m~'^ exp{pm) = {m*)~'^^{m/m*). 
Here, the scaling function fl>(x) = x~'^ exp(—x) with 
m* ~ p~^ and power law exponent r = 2 -|- l/i5 with 
1 < T < 2 (as (5 < —1), leading to a relation 5{t — 2'] = 1. 
The scaling form was numerically observed in [34| . In¬ 
terestingly, the borderline case with (5 = — 1 generates 
gamma distributions, which are found in a broad class 
of mass transport processes 3l| and have been also ob¬ 
served in a limiting case of conserved-mass aggregation 
models studied in Ref. M. 


B. Nonconserved Mass Aggregation Models 

In this section, we discuss a nonconserved version of 
the mass aggregation models where systems can exchange 
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mass, though weakly, with environment. In this case, in 
addition to the earlier defined two processes (1) diffusion 
and (2) fragmentation of masses, a particle now can be 
adsorbed with rate q and desorbed at a site with rate p, 
provided the site is occupied, where p, g —?► 0 (i.e., weak 
exchange) with the ratio r = q/p finite. Due to adsorp¬ 
tion and desorption processes, total mass in the system 
is not conserved. This model is related to several mod¬ 
els studied in the past for finite p and q 
Interestingly, in the limit oi p,q —^ 0, mass fluctuation 
in a nonconserved model can be obtained from the occu¬ 
pation probability of a site in its conserved version (i.e., 
p = q = 0) 1^. Let us first define, in the space of 

total mass M, a transition rate Tm+i,m from mass M to 
M + 1. In the steady state, the probability P{M) that 
the system has mass M satisfies a balance condition 


P{M)Tm+i,m = P{M + 1)Tm,m+i 
where the mass distribution P{M) can be obtained as 


P{M + l) 


11^T(M + I^M) 


P(0). 


As the ratio of transition rates can be written as 


Tm+i,m _ q 
Tm,m+i pS{p) 


(30) 


where S{p) is the occupation probability and p = M/V, 
the distribution function can be written, upto a normal¬ 
ization factor, as 


P{M) OC ^ g-V /J” dp[/i(p)-/io] 


(31) 


where po = lii(9/p) an effective chemical potential and 
f{p) = f dpp{p) = f dp In S(p) an effective free energy 
(canonical) density function. The steady state mass 
density as a function of adsorption to desorption ratio 
r = qjp can be obtained by minimizing the grand poten¬ 
tial or the large deviation function for the density fluctu¬ 
ation /i(p) = /(p) — PoPj leading to the relation 5(p) = r 
(for details, see Appendix C). 

Till now, the analysis is exact. However, it may not 
always be possible to exactly calculate the occupation 
probability 5(p). For the purpose of demonstration, let 
us proceed by considering a model with diffusion and 
fragmentation rate as in Case I. We obtain an approx¬ 
imate expression, obtained within mean-field theory, of 
the occupation probability (see Appendix C) 


S[p) = 


P(1 - P) 
(1 + p) ■ 


Then, Eq. [21] implies the subsystem mass distribu¬ 
tion having a form Py{m) oc w„(m) exp(p,m) and con¬ 
sequently a FR relation as in Eq. [3] follows. Then, for 
p < Pc or equivalently for r < Tc, one can immediately 
calculate the scaled variance as 


<(p) = (^ 


-1 


p(I-p)(I + p) 
(I-2p-p2) ’ 


(32) 


where critical density Pc = ^/2— 1. The variance in non¬ 
conserved case is different from that in the conserved- 
mass case, implying that the canonical a nd g rand canon¬ 
ical ensembles are not equivalent [s^, l38l | However, 
the nature of singularity in the variance remains the 
same near criticality where cr^(p) (pc ~ p)~" with 
n = 1. Therefore the additivity property leads to the 
same power law scaling in the single-site mass distribu¬ 
tion Pi (to) ^ m~'^ exp{pm), for large to, where t = 5j2 
and p, = po — ln5(pc) = ln(r/rc) with Tc = 5(pc). 

The above results are consistent with what was found, 
on the mean-field level, for general p and q in ‘in-out’ 
model dH - a special case of the above nonconserved 
model with w = 0. One can interpret the results in the 
light of equilibrium BEG: The critical density signifies 
that, for r > Tc = S{pc) = 3 — 2\/2, there is a conden¬ 
sate as in the BEG. In the grand-canonical setting (i.e., 
with no mass conservation), that would imply a phase 
with a diverging mass density, similar to the ‘Takayasu 
phase’ where mass density actually diverges. For p and q 
finite, form of the subsystem mass distribution as written 
in Eq. [31] remains the same, but only that the expres¬ 
sion of 5(p), due to the presence of spatial correlations, 
is different. However, the similarity with the BEG still 
persists. 


IV. SUMMARY AND CONCLUDING 
PERSPECTIVE 

In this paper, we argue that an additivity property can 
possibly explain why simple power-law scaling appears 
generically in nonequilibriuin steady states with short- 
ranged correlations. We demonstrate that the existence 
of a fluctuation-response relation, a direct consequence 
of additivity, with a singular response function leads to 
power-law distributions with nontrivial exponents. The 
simplest form of the singularity, a simple pole, gives rise 
to the exponent 5/2, which was often observed in the 
past in apparently unrelated systems. We substantiate 
the claims by analytically calculating the response func¬ 
tion, which diverges as critical point is approached, in 
paradigmatic nonequilibrium mass aggregation models 
and the corresponding single-site as well as subsystem 
mass distributions. Most remarkably, the analysis, being 
independent of dynamical rules in a particular system, 
equally extends to critical properties in equilibrium and 
nonequilibrium. 

Thermodynamic characterization of phase coexistence 
in driven systems is a fundamental problem in statistical 
physics and has been addressed in the past [2l[2i,|3i- 
l43l |. either numerically or analytically only for exactly 
known steady-states mostly having a product measure. 
From that perspective, it is quite encouraging that, even 
when steady-state weights are a priori not known, our 
analytical method not only gives insights into the steady- 
state structure but can also be applied to identify a chem¬ 
ical potential, which equalizes in the coexisting phases 









and vanishing of which at the criticality gives rise to pure 
power laws. 

Note that, in our formulation, the mass distribution 
functions, though approximate, have been calculated 
solely from the knowledge of the variance. This formula¬ 
tion is perhaps not surprising in equilibrium where free 
energy function (or entropy, for an isolated system) es¬ 
sentially determines fluctuation properties of a system. 
However, in nonequilibriuin scenario, it is a-priori not 
clear that such equilibrium thermodynamic approach can 
indeed be applied in systems having a steady state with 
nontrivial spatial structure. Here, it is worth mention¬ 
ing that one requires, in principle, all the moments to 
specify a probability distribution function. However, ad¬ 
ditivity property, provided it holds, puts a strong con¬ 
straint on the mass distribution function Py (m) through 
a fluctuation-response relation and thus helps to uniquely 
determine Py(m), only from the knowledge of the vari¬ 
ance as a function of density. 

We believe that our analysis, being based on a general 
thermodynamic principle, would be applicable in many 
other driven systems where phase coexistence is known 
to occur (e.g., in active matters Qli^). As a concluding 
remark, we mention that additivity property is expected 
to be quite generic for systems having short-ranged cor¬ 
relations and, therefore, it would be indeed interesting to 
actually verify additivity, through the predictions con¬ 
cerning fluctuations, on a case-by-case basis. Also, it 
remains to be seen whether the principle of additivity 
can be extended to systems having long-ranged spatial 
correlations, at least in the cases where the strength of 
these correlations are weak. 
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APPENDIX 


APPENDIX A: CALCULATION OF VARIANCE 
IN CONSERVED MASS AGGREGATION 
MODELS (CMAM) 

We define here a class of conserved mass aggregation 
models (CMAM) on a one dimensional lattice with pe¬ 
riodic boundary and calculate the variance of mass at a 
single site in the steady state, assuming that the additiv¬ 
ity property (Eq. 1) holds. For, simplicity, we consider 
only the discrete-mass cases. 

The mass at each site undergoes either diffusion (where 
whole of the mass is transferred to either of neighbouring 
sites) or chipping, with certain transition rates; in the 
models considered below, there are two types of chipping 
process. The diffusing mass or the chipped-off mass are 
coalesced with the mass at either of the neighbouring 
sites with a pre-assigned rates. In this process, the total 
mass of the system is conserved. 

Provided a site i is occupied, particles hop to either of 
the nearest neighbour sites according to the transition 
rates specified below. 

A. Diffusion with rate 1: All particles at a site i 

hop with rate 1 to left or right, i.e., > 0 and 

mi±i -)■ TOiii -I- m*. 

B. Chipping with rate wi: This chipping process 
involves a particle at site i being chipped off and thrown 
to left or right neighbour, i.e., rrii {irii — 1) and 
TOiii -)■ TOiii -I- I. 

C. Chipping with rate W 2 : This chipping process 
involves mi — 1 particles going to either left or right 
neighbour and the rest of the particles remaining at site 
f, i.e., mi ^ I and mi±i mi±i + mi — 1. 

Assuming transition rates are Poissonian, we have 
the following stochastic update rules where mass 
mi(t + dt) at time t + dt takes a particular value, depend¬ 
ing on mass mi[f) at time t, with certain probabilities 
as shown below. 

Loss terms at site i: 


value: probability: 

0 dt 

mi{t) - I -h widt 

1 W2dt. 


In this Appendix, we provide the details of the calcu¬ 
lations to obtain the mass distributions, using additivity 
property, in mass aggregation models (both conserved 
and nonconserved versions) which were studied over the 
last couple of decades. The generalized models intro¬ 
duced here cover some of those studied in the past and 
their variants [il,[il,[iiill,[2i,|33 


Gain terms from {i — 1)*^ site: 


mi(t + dt) 


value: prob.: 

m^{t) + m^-l{t) 

mi(t) -t- I T ’ 

m^(t)+m^-i(t) - 1 + W 2 f. 
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Gain terms from {i + 1)*^ site: 


mi(t + dt) 


value: prob.: 

m,{t) + 

mi{t) + 1 - wif 

mi{t)+mi+i{t)-l + 6mi+i(t),o W2^. 


Now we define the occupation probability 
((1 — Smj,o)) = S{p), i.e., the probability that a 
site is occupied. We deal with steady-state averages 
throughout. We assume that the additivity property 
(as in Eq. 1) is valid at single site level and therefore 
n-point (n > 2) correlation factorizes. 


Mass remains unchanged at site i: 

value: prob.: 


mi{t + dt) = 


i(t) (1 — 2dt — 2widt — 2 w 2 dt). 


n-th moment equation: The time evolution of 
n-th moment (to”) can be written as 


(TO”(t -f dt)) = {mi{t)) = {[mi{t) -1 + Sm,{t),or)widt + {[m^{t) + TOi_i(t)]”)y -f ([1 - 6 mi{t),or)w 2 dt 

+{[m^(t) + 1 - + {[m^it) + mi-i(t) -1 + (5„,_i(t).o]"')w'2y + ([TOi(t) -f TO*+i(t)]”)y 

+ {[mi{t) -I- 1 — + ([TOi(t) -I- TOi+i(t) -1-1- <^mi+i(t).o]"')^2y 

-|-(TO”(t))(l — 2dt — 2widt — 2w2dt), (33) 


which, in the steady state where {rn^{t + dt)) = {m2{t)), 
gives a BBGKY hierarchy where n-point correlations 
are coupled to (n -|- 1 (-point correlations. To get a 
closed set of equations for the moments, we use the 
factorization property of n-point correlations. As 
mentioned in the paper, the mass distributions are 
solely obtained from the response function (or the 
variance of the mass distribution) and therefore we are 
interested in only calculating the variance, or equiv¬ 
alently the second moment, which can be done as follows. 


2"*^ moment equation: If we put n = 2 in the 
above equation, the second moment {mf) however 
cancels out from the above equation. Using factorization 
of two-point correlation, i.e., (rrnmj) ~ for i ^ j, we 
get an expression for the occupation probability S{p) as 
a function of mass density p, 

p^{l +W 2 ) = w+{p-S) - w-pS, (34) 


where w± = w\±W 2 - This gives 

S{p) = (1 +^2)p^ 

W+ W-P 


(35) 


3'’'^ moment equation: Similarly, for n = 3, we 
obtain an equation where the third moment (to?) 
cancels out and we actually get, using factorization of 
both two-point and three-point correlation, a relation 
for the second moment 


rc+(l -I- 5) — 2w2P 
u>+ — 2(1 -I- W 2 )p — W-S 


Using the expression of occupation probability in Eq. 1351 
we obtain 

2 _ w\ + 2w+w-p - (w+ + 3 wiW 2 - wl)p^ 

^ ^ w\- 2 w+{l+W 2 )p-'W-{l+W 2 )p‘^ 

which leads to the desired expression of the variance as 
a function of density, 

2 , ^ _ w\p + w+{wi - 3 w2)p^ _ 

^ w\-2w+{l + W2)p- W-{l+W2)p^ 

(w+ - W 1 W 2 + 3wj)p^ -K W-(l -H W 2 )p '^ , . 

— 2w+{l + W 2 )p — W-(l -I- W 2 )p^ 

The variance a^ip) has a singularity at p = pc, i.e., it di¬ 
verges at a critical density p = pc, which can be obtained 
by putting the denominator of Eq. [3H] zero and solving 

— 2w+{l + W 2 )pc — w-(l -I- W 2 )p^ = 0. (39) 

This gives a simple pole at the critical density 


Pc 


w+ 

W- 



W- 

1 -I- W2 



(40) 


Nonequilibrium free energy function can be calculated 
by integrating nonequilibrium chemical potential w.r.t. 
density p, 

dip) = ^ ^ = j P^P)<^P- ( 41 ) 

The function A„(s) = — ln?i;(s), which is the Legendre 
transform of the free energy density /(p), can be obtained 
as given below. 


(to^) = p 


(36) 


\{s) = v[f{p*) -H sp*]. 


(42) 
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where p* is the solution of 


s = -p{p*). 


APPENDIX C: CALCULATION OF MASS 
DISTRIBUTION IN THE ABSENCE OF MASS 
(43) CONSERVATION 


APPENDIX B: CALCULATION OF MASS 
DISTRIBUTION IN THE CONSERVED MASS 
AGGREGATION MODEL 


As shown in the paper, the probability distribution 
function P{M) of total mass M can be written, up to a 
normalization factor, as 


Here we provide the essential steps of the calculations 
to obtain single-site (i.e., w = 1) mass distribution func¬ 
tion Pi{m) oc tyi(m) exp[/x(p)m] where wi{m) is the 
single-site weight factor and p,{p) is a chemical potential. 
We first analyse the behaviour of Ai(s) near the singu¬ 
larity at s = Sc by expanding p{p) and /(p) near critical 
density in the power series oi p — Pc where p — Pc < 0 is 
small, 

m(p) = KPc) + - Pcf + ■■■ (44) 

fip) = f{Pc) + p{Pc){p-Pc) + Pc)^ + ■ ■ ■ 

where we have used Eq. [41] and p'{pc) = f"{Pc) = 0 
(see Fig. |T|). Using Eq. |33| in Eq. 011 and substituting 
s -f p(pc) ~ -p"(pc)(p* - Pc)^/2, we get 


where Sc = —p(pc) and p"{pc) < 0. Therefore Ai(s) = 
/(p*) -I- sp* near s = Sc, in the leading order of (s — Sc), 
can be approximated as 


Ai(s) 


f{Pc) - Sc{p* - Pc) + ^-^^{p* - Pc)^ 

Al(Sc) -I- p*(s — Sc) -I- ^ (p* — Pc)^ 

Oo -I- ai(s - Sc) -I- a 2 (s - Sc)^'^^ 


I ^ 

-I- sp 

(46) 


where oq = Ai(sc) = /(pc) + ScPc, oi = Pc and 02 = 
— (2/3)a/2/|p"(pc)|. The inverse Laplace transform of 
the weight factor wi{m) can be written as 


u)i(s) = ~ e-“«[l-ai(s-Sc)-a2(s-Sc)3/2] (47) 

which, for m ^ 1, translates into 

gScm 

wi{m) -(48) 

vnP' 

Consequently the mass distribution can be written as 


Pi(m) 


^Scm 


, , p-(a+iJ.o{pc))m 

pP{p)m _ ^ _ ^{po(P)+a)m 

^5/2 

Pi{m) - (50) 

nnpl ^ 


Note that effective chemical potential /i(p) = po(p) — 
Po{pc) is zero at the critical density pc = (v^ — 1). The 
mass distribution in Eq. (SOj is precisely what was found 
in [T^ at p = Pc and describes the simulation data re¬ 
markably well (see Fig. 1). 


P{M) = const. X e"^ -fo (51) 

Now, if we assume that the joint mass distribution 
P[{mi}] has a product form on single-site level {v = 1), 
i.e., product of single-site mass distribution function 
Pin^i), 


V 

'P[{mi}] = qp(mi), (52) 

i=l 

the probability distribution function P{M) of mass M in 
the system can be written as 

p{M)= n 

i=l 

From the Laplace transform U(s) = 

f dMP{M) exp(—sM) = [p(s)]'^ of the mass dis¬ 

tribution P{M), the Laplace transform p{s) = 
J dmip(mi) exp(—smi) of single-site mass distribu¬ 
tion p{m) can be written as 

p(s) = const. X (54) 

where 

Ai(s) = infp[/i(p) -I- sp]. (55) 

Here we have used inverse transform 

P(s) = const. X J (56) 


dmipijrii) 


5\M- 


(53) 


which has been obtained from Eq. 01] and where grand 
potential or the large deviation function for density fluc¬ 
tuation h{p) = /(p) —pop = /p^[p(p) —pojdp and chemical 
potential p(p) = ln5(p) = ln[p(l —p)/(l-|-p)], as given in 
the paper. Note that the function 5(p) is the occupation 
probability in the conserved mass aggregation model and 
has been obtained by putting ici = 1 and iU 2 = 0 in Eq. 

ISSl 

Now the function Ai(s), Legendre transform of grand 
potential /i(p), can be written as 

Xi{s) = h{p*) + sp\ (57) 


where p* is the root of the equation d[h{p) sp]/dp = 0 
or p(p*) — po + s = 0, i.e., p* is the root of 


In 


P*{l-P*) 

l+P* 


= po - s. 


(58) 
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The critical density is obtained by putting scaled variance 
as cr^{p) = {dp,o/dp)~^ = oo or l/a'^{p) = 0, 


(1 - 2pc - pI) 

Pc(l - Pc){l+ Pc) 


(59) 


to obtain 


Ai(s) 


(S - Sc) 

ao + ai(s — 


W'iPc)\ 

2 


iP* 


Pc)" 


Sc) + 02(3 - 


(61) 

(62) 


and thus pc = \/2— 1. In the macrostate (most probable 
state), we have S{p) = r, implying that the critical den¬ 
sity is related to the ratio r = q/p through 5(pc) = ?’c- To 
obtain the large-mass behaviour, we expand p{p) around 
P ~ Pcj 

dip) = P-iPc) + ^ - Pc)^ (60) 


in leading order in {p* — pc) where Sc = po—p{pc), leading 
to the desired result in the paper, 

Pirn) ^ 1 (63) 
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